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ABSTRACT 

We present a general scheme for constructing Monte Carlo realizations of equilibrium, col- 
lisionless galaxy models with known distribution function (DF) /q. Our method uses impor- 
tance sampling to find the sampling DF fg that minimizes the mean-square formal errors in 
a given set of projections of the DF /q. The result is a multi-mass A^-body realization of the 
galaxy model in which "interesting" regions of phase-space are densely populated by lots of 
low-mass particles, increasing the effective N there, and less interesting regions by fewer, 
higher-mass particles. 

As a simple application, we consider the case of minimizing the shot noise in estimates 
of the acceleration field for an A^-body model of a spherical Hernquist model. Models con- 
structed using our scheme easily yield a factor ~ 100 reduction in the variance in the central 
acceleration field when compared to a traditional equal-mass model with the same number 
of particles. When evolving both models with a real A^-body code, the diffusion coefficients 
in our model are reduced by a similar factor. Therefore, for certain types of problems, our 
scheme is a practical method for reducing the two-body relaxation effects, thereby bringing 
the A^-body simulations closer to the collisionless ideal. 
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1 INTRODUCTION 

There are two types of A^-body simulations in stellar dynamics. In 
collisional simulations each of the A^ particles represents an indi- 
vidual star. This type of simulation is most often used to model 
the evolution of star clusters in which discreteness effects, such as 
two-body relaxation, are important. 

When it comes to modelling galaxies, however, the number 
of stars is large enough and the dynamical time is long enough 
that these discreteness effects are usually unimportant. In the limit 
of a very large number N of bodies, stars and dark matter par- 
ticles move in a smooth mean-field potential <^{x\t) and behave 
as a collisionless fluid in six-dimensional phase-space (Binney & 
Tremaine 1987, BT87), the (mass) density at any point (x, v) be- 
ing given by the distribution function (hereafter DF) fix, v;t). The 
time-evolution of the DF is described by the Collisionless Boltz- 
mann Equation (hereafter CBE). Therefore, in a collisionless N- 
body simulation the TV particles do not correspond to real stars; 
instead they provide a Monte Carlo realization of the smooth un- 
derlying DF, from which one can estimate the potential ^{x,t).'By 
integrating the orbits of these particles, one is solving the CBE by 
the method of characteristics (Hernquist & Ostriker 1992, H092; 
Leeuwin Combes & Binney 1993, LCB1993) 

In reality, no simulation is perfectly collisionless because 
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Poisson noise in the estimates of ${x,t) inevitably leads to nu- 
merical diffusion in particles' orbits. To reduce this noise, it is im- 
portant to make N, the number of particles in the simulation, as 
large as possible. Unfortunately, the cost of running an A'-body 
code scales at least linearly with A^, so increasing A'' also makes 
the simulation more expensive to run. The good news is that al- 
ternative, more sophisticated weapons are available for use in the 
fight against small A^ limitations. A collisionless A^-body code is 
essentially a Monte Carlo method and so should be amenable to 
well-known variance-reduction methods such as importance sam- 
pling (e.g.. Press et al 1992). 

In this paper we present a generally-applicable, essentially 
model-independent method for constructing A^-body realizations 
of isolated model galaxies in equilibrium, suitable for use as initial 
conditions (hereafter ICs) in collisionless simulations. Our scheme 
uses importance sampling to find a sampling DF fs that minimizes 
the mean-square uncertainty in a chosen set of projections of the 
DF fo . For example if modelling bar evolution, one might be most 
interested in following the detailed evolution of the DF around the 
strongest resonances. It is natural then to try to increase the sam- 
pling density near these regions by populating them with lots of 
low-mass particles. Outside these interesting regions, however, one 
must also have enough particles to maintain accurate estimations 
of the force field which governs the evolution of the system as a 
whole. 

The paper is organized as follows. After reviewing the basics 
of Monte Carlo integration and the connection between CBE and 
A^-body simulations, we explain our multi-mass formulation in sec- 
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tion 2.4. With the notable exception of some heuristic multi-mass 
schemes (e.g., Sigurdsson, Hemquist & Quinlan (1995), hereafter 
SHQ95, Weinberg & Katz (2007), Sellwood (2008) and Zemp et al. 
(2007)), most other IC-generation schemes have used equal-mass 
particles. In section 3 we give an example of using our scheme to 
suppress fluctuations in the monopole component of acceleration 
in a spherical galaxy model. We calculate formal estimates of the 
noise in A^'-body models constructed the equal-mass scheme (sec- 
tion 3.1.1) and SHQ95's method (section 3.1.2) and compare them 
to our our own scheme in section 3.1.3. In section 3.2 we test how 
well our realizations behave in practice when evolved using a real 
A'^-body code. Finally, section 4 contains a summary of the prereq- 
uisites for our scheme, along with possible scientific applications. 



2 FORMULATION 

2.1 Monte Carlo integration 

For later reference we recall some of the basic ideas (e.g.. Press et 
al 1992) in using Monte Carlo methods to evaluate integrals, such 

as 



fdV, 



(1) 



of a known function / over a domain D. We first consider the case 
where D has unit volume: dV — 1. Then, given points, 
xi . . .xn, drawn uniformly from D we can estimate 
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as A*' ^ cxD. 

Now let us relax the assumption that D has unit volume and, 
instead of drawing points uniformly from D, let us take A'^ points 
drawn from a sampling distribution fs. We assume that fs is nor- 
malised: fs AV = 1. Making a straightforward change of vari- 
ables and using the result above, it follows that the integral (1) can 
be estimated as 



1 ^ 
N ^ 



and that the variance in this estimate is approximately 



dV 



(4) 



(5) 



2.2 A^-body simulations and the CBE 

The following description of the connection between A'-body sim- 
ulations and the CBE borrows heavily from LCB93. We assume 
that the galaxy has total mass M, = 1 and that the mass density of 
stars in phase-space is given by a DF f{x,v;t) = f{w;t), where 
w = {x,v), normalised so that 

f{w)d^w^l. (6) 

The evolution of the DF is governed by the CBE, 

^-0 (7) 
dt dx dx dv 



It conserves phase-space density, so that 

f{w{t);t)^f(wo;0), (8) 

where w (t) is the path traced by an individual particle, with wo = 
w{t = 0). As H092 and LCB93 point out, in a colhsionless A^- 
body simulation one is solving the CBE for these w{t) by integrat- 
ing the characteristic equation, 

dt da; dv 
1 V a 

and using Monte Carlo integration to estimate the acceleration 

^('^'''^d^^,'. (10) 



x'\ 



From (4) it follows that 



a{x;t) ~ -GV^ ■ 



(11) 



corresponding to a distribution of A'^ point particles with masses 

N fs(wi;t) 

These m; clearly depend on the choice of sampling DF fs. The 
simplest choice is fsiw; t) — f{w;t) — f{wo; 0), in which case 
all particles have equal masses rrii — 1/N. However, one is free to 
tailor the choice of fs to suit the particular problem under study. 

The singularities in (11) at a; = Xi yield estimates of a{x) 
that suffer from unacceptably large scatter; in fact, they corre- 
spond to the direct accelerations appropriate for a coUisional A'- 
body code! So, in practice colhsionless simulations do not use (11) 
directly, but instead obtain a{x) using techniques (e.g., softened 
force kernels, grid methods or truncated basis-function expansions) 
that reduce the scatter by removing the singularities. More gener- 
ally, eq. (1 1) provides only the most simple-minded estimate of the 
integral (10), and one has some leeway in how one reconstructs 
f{w;t) from the discrete realization furnished by the A'' parti- 
cles. Of course, the reliability of any sensible reconstruction will 
be wholly dependent on how well the DF is sampled. 

2.3 Observables 

What constitutes a "good" choice of sampling density /s? The DF 
/ is a high-dimensional probability density and itself is not measur- 
able. We are usually only interested in coarse-grained projections 
of the DF 



f{w)Q,{w)d''w, 



(13) 



where the kernels Qi{w) are some functions of the phase-space co- 
ordinates {x,v). For the purposes of the present paper, we consider 
a "good" sampling scheme to be one that minimizes the uncertainty 
in the estimates of some given set of (Qi). Apart from some general 
guidance, we do not address the question of how best to choose 
these Qi, which usually requires some experience of the particular 
problem at hand. 

We now give some examples. It is helpful to introduce the 
indicator function 



lv{w) 



if 10 e 1/ 

otherwise. 



(14) 



Then a particularly simple but important choice of kernel is 
Q,{w)^lvAw), (15) 
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for which (Qi) measures the mass inside a volume Vi. For many 
problems one might choose some of the Vi to surround important 
resonances in phase-space, so that (Qi) measures the phase-space 
density around the resonances. With appropriate choices of projec- 
tion kernel Qi, the expression (13) includes quantities such as the 
galaxy's density profile, its velocity moments or even its projected 
line-of-sight velocity distributions. 

More fundamentally, an A'^-body model should provide a good 
estimate of the galaxy's acceleration field. Therefore we recom- 
mend that many of the (Qi) be used to measure at least the 
monopole component of the galaxy's acceleration field at a range of 
points. This can be achieved using (15) with spherical volumes Vi 
centred on x — for a range of radii ri, encompassing all velocity 
space for |x| < ri. Similarly, one can include higher-order mul- 
tipole components of the galaxy's acceleration field by choosing 
a slightly more complicated projection kernel Qi (see equation 35 
below). 



We partition phase space into n/ cells and write Tj ^ for the 
phase-space volume enclosed by the j^^ cell. We parametrize fs as 



1., 



(21) 



3 = 1 



so that within the j' phase-space cell fa is given by fo{w)/aj. 
For the equilibrium models considered, it is natural to choose tj to 
be cells in integral space. Substituting this fs into (19) yields 



'^ajHi 



where 



If we further define 



(22) 



(23) 



2.4 Optimal sampling scheme 

The problem we address in this paper is the following. We wish 
to construct an equilibrium A*'-body realization of a galaxy model 
with some known DF /q. Specifically, we seek ICs that faithfully 
represent some projections. 



foQi d w, 



(16) 



of this DF, for a set of ng kernels Qi{w). What is the "best" choice 
of sampling DF fs given this fo and choice of kernels Q, ? 

More formally, from (5) the uncertainty in a Monte Carlo esti- 
mate of (Qi) obtained using A'^ particles drawn from the sampling 
distribution fs is given by 



Var{Q,) 



TV 



fs' 



) ,, d W 



(17) 



Notice that, unlike most introductory textbook examples of Monte 
Carlo methods, we have ?iq such estimates but just one fs. We 
seek a normalised sampling DF fs that minimizes the mean-square 
fractional uncertainty 



(18) 



where SQi, the formal relative uncertainty in a measurement of Qi, 
is given by 



(SQ, 



2 _ Var (Q.) 



(19) 



Of course there are many other possible measures of the "good- 
ness" of some choice of fs . 

One can immediately use the Euler-Lagrange equation to 
show that choosing 



fs(w) CC fo{w)}_^ ' 

[Qi) 



(20) 



extremizes (18), the proportionality constant being set by the con- 
straint that fs should be normalized, J fs d^w = 1. This direct so- 
lution is flawed, however, since for most interesting choices of Qi 
the resulting fs depends on orbit phase; using this fs the masses of 
particles sampling a given orbit would vary along the orbit! There- 
fore in practice we use a slightly less direct approach. 



then the mean-square fractional uncertainty (18) becomes 



ng 



j=i 



(24) 



(25) 



Our goal is to find the coefficients aj that minimize this S, subject 
to the constraint that the resulting fs be normalised. The normali- 
sation constraint is that 



where 



/() d'^w. 



(26) 



(27) 



Using the method of Lagrange multipliers, the coefficients of the 
"best" sampling DF obtained by minimizing (25) subject to the 
constraint (26) are simply 




(28) 



which is just the direct solution (20) in disguise, but averaged over 
the phase-space cells Tj and correctly normalized. This averaging 
means that the resulting fs will be smooth, provided that none of 
the kernels Qi pick out specific regions of integral space. 

Substituting the fs given by (21) into (12), we have that par- 
ticles in phase-space cell r, have masses rrij — aj/N. One can 
therefore easily impose additional, direct constraints on the masses 
of particles within a subset of the phase-space cells r, ; simply re- 
peat the minimisation of (25) subject to (26) while holding the rele- 
vant subset of the aj fixed at the desired values. For example, when 
generating an A'^-body realization of a dark-matter halo model in- 
side which one intends to embed a disk of light particles, one might 
want to ensure that those halo particles passing through the disk 
have the same mass as the disk particles. Of course, a more pedes- 
trian approach would be to introduce additional kernels Qi to pick 

^ Note that we use V to denote subvolumes of phase-space be used in the 
calculation of the projections (13) of the DF fo, and r for the subvolumes 
used in the discretization of the sampling DF fs . 
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out the relevant parts of integral space. We caution, however, that 10' 

we have not tested how well such a "bumpy" fs would work in 

practice; the tests we present later all involve smoothly varying 

sampling distributions. ^^o 



2.5 ICs for iV-body model 

Together with /o, the coefficients aj completely determine the sam- 
pling DF fs of the form (21). In particular, it reduces to the con- 
ventional equal-mass case when all aj = 1. 

We apply the following sequence of steps A'^ times to draw 
particles from this fs, thereby constructing an A'^-body realization 
of the galaxy model: 

(1) Choose one of the Uf cells at random, the probability of 
choosing the j^^ cell being given by Ij/aj. Let i be the index of 
the chosen cell. 

(2) Assign a mass m, = fo{wi) /N fs{wi) = ai/N to the 
particle. 

(3) Within the i^^ cell, draw a;, from its density distribution. 
Pi = J folrid^w. For the special case of a spherical galaxy, one 
can precompute the cumulative mass distribution Mi{r) for each 
of the 71/ cells and use this to draw a radius r^, followed by angles 
9i and cjii . 

(4) Use an acceptance-rejection method to draw Vi from 
fo{xi, v) at this fixed value of Xi. 



10-^ ; 

10"" 10"^ 10"^ 10"' 10° 10' lo'^ 
r/a 

Figure 1. Formal relative errors 5M = (Var (M))^/^/ (A/) (eq. 19) in 
the monopole component of the potential of a Hemquist model (eq. 29) 
constructed using the same number A'^ = 10^ of particles, but drawn from 
different sampling DFs. The heavy solid curve plots results for our tailored 
sampling DF. For comparison, we also show results for the conventional 
equal-mass scheme (light solid curve) and Sigurdsson et al.'s (1995) heuris- 
tic multi-mass scheme (dashed curve). 



3 AN EXAMPLE 

In this section we use a simple galaxy model to demonstrate our 
scheme. Our galaxy model is spherical and isotropic, with density 
profile (Hemquist 1990) 



p{r) 



M, 



2nr{r + a)^ 



(29) 



total mass Ah and scale radius a. By Jean's theorem, the model's 
DF f{x, v) depends on {x, v) only through the binding energy per 
unit mass £. Hernquist (1990) gives an expression for /(£). 

We want to construct an A^-body realization of this model 
that minimizes the mean-square error in the monopole component 
of the acceleration averaged over many decades in radius, from 
^min ~ 10~*a up to Tmax = lO^a. To achievc this we choose 
kernels Qi — Iv; (r) that measure the mass enclosed within a se- 
quence of 25 spheres centred on the origin, with radii ri spaced 
logarithmically between r-min and rmax- We use (22) to calculate 
the formal uncertainty SMi in the enclosed mass for a range of dis- 
cretized sampling densities of the form (21), including (28). 

To implement this, we first of all partition integral space 
(£,j'^) onto a regular 71/ — ne x nx grid. The ng energies £j are 
chosen to match the potential £j = ^{rj) with Vj logarithmically 
spaced between lO^^a and lO'^a. At each £j, there are nx values 
of Xjk running linearly from to 1, where Xjk = Jk{£j)/ Jc{£j) 
is the orbital angular momentum normalised by the circular angu- 
lar momentum at energy £j ^. These choices ensure that our fs 
samples well the interesting parts of phase-space. For the calcula- 
tions below we take ne x nx = 200 x 100, although a coarser 
grid (e.g., 50 x 25) would suffice. Having defined our projection 
kernels Qi — Ivi, we use (16) to calculate the expected values of 



^ X = J(£)/ Jc{£) is a measure of an orbit's circularity; orbits with 
X = I are perfectly circular, while those with X = are perfectly radial. 



enclosed mass {Mi) and the ancillary quantities Hij (from eq. 23), 
and use these to obtain the formal uncertainties 5Mi in (22). 



3.1 Comparison with other schemes: formal errors 

Before applying our method, we study two other schemes: the con- 
ventional equal-mass scheme and the multi-mass scheme of Sig- 
urdsson, Hernquist & Quinlan (1995). 



3.1.1 The conventional equal-mass scheme 

The most common (albeit implicit) choice of sampling density is 
fs = /o, which corresponds to setting all Oj = 1 in our equa- 
tion (21). All particles then have the same mass. For our ex- 
ample Hemquist model the fraction of particles within radius r is 
r^/(a+r)'^, so that less than 1% of the particles are within 0.1a. As 
shown in figure I, for this fs the formal uncertainty 5M{r) rises 
steeply towards the centre. Although this scheme produces accu- 
rate estimates of the galaxy's potential outside the scale radius a, it 
performs poorly in the interesting r^^ central density cusp. 



3.1.2 Sigurdsson et al.'s multi-mass scheme 

Sigurdsson, Hernquist & Quinlan (1995) have used an interesting 
heuristic scheme to improve the resolution of A^-body models near 
galaxy centres. In effect, they use an anisotropic sampling function 
of the form (21) with coefficients 



a(T) = B X 



if T-peri < a, 

otherwise. 



(30) 



where rpcri(T) is the smallest pericentre radius of any orbit from 
the phase-space cell r, and the constant B is chosen to normal- 
ize fs. When the parameter A = 0, then aj — 1 and the sampling 
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Figure 2. Contour map of logjg (/s/Jo) for the optimal muli-mass sam- 
pling scheme of section 3.1.3; notice the strong enhancement in low angular' 
momentum and large energy (small rg ) region. These correspond to orbits 
with small peri-centre radii. 



DF fs is identical to /q. Increasing A improves the sampling of 
the cusp by increasing the number density of particles having peri- 
centres Tperi < a. Consequently, as r — > the number density of 
particles rises more rapidly than the mass density, permitting better 
resolution in the centre. To balance this increase in number den- 
sity, each particle is assigned a mass fo/Nfs = aj/N so that the 
phase-space mass density is still given by the desired fo . 

The dashed curve in figure 1 shows the formal error SAI(r) in 
our implementation of their scheme for A = 1. Their scheme does 
much better than the conventional equal-mass scheme at small radii 
r ^ a, at the cost of a slightly noisier monopole at r > a. 



3.1.3 Our scheme 

It is encouraging to see that SHQ95's multi-mass scheme does, to 
some extent, improve mass resolution at small radius. However, as 
shown in figure 1, 5M at r = 10~''a is still almost two orders of 
magnitude larger than at r = a. Can we achieve even better results 
by carefully designing an fs that generates a flat 5M{r) across a 
large range of radii? 

The fs given by our optimal choice of coefficients (28) is 
plotted in figure 2. It is qualitatively similar to SHQ95's results, in 
the fact that it samples densely the low-angular momentum parts of 
phase-space. The detailed shape of the function is different, how- 
ever, and the thick solid curve in figure 1 shows that our scheme 
provides much better estimates of the monopole components of the 
acceleration at small radii; in fact, the formal error 5M{r) varies 
by only a factor ~ 4 over six decades in radius. 



3.2 TV-body realizations 

Figure 3 shows the spectrum of masses obtained using the algo- 
rithm detailed in §2.5 to draw — 10® particles from this optimal 
fs. Unlike the conventional scheme which would give all particles 
the same 10* mass if assuming A/* = 10^" Mq, our multi- 
mass scheme assigns a range of masses between 10~^A/o and 
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Figure 3. Histogram of particle masses from an A'^ = 10*^ multi-mass real- 
ization of a Hernquist galaxy, scaled to a total mass A/* = 10^'' Mq. The 
span of 8 decades in mass gives sub-solar mass resolution in "interesting" 
regions of phase space. In contrast, in an equal-mass realization all particles 
would have mass lO'^M© (thick solid line). 

lO^Afe (8 decades), with many low-mass particles in the central 
region. 

As a simple sanity check of our formal estimates of the errors 
in the monopole, we count the mass of particles within the same 
spheres Vi used to calculate 5 Mi. The deviations from the mass 
profile of the target Hernquist model are consistent (figure 4) with 
the expected values of 5M from equation (22). 

Ultimately, the purpose of our sampling scheme is to im- 
prove the numerical modelling of collisionless galaxies close to 
equilibrium using full A'^-body integrations. To test how well our 
scheme succeeds at this task, we use the particle-multiple-mesh 
code Grommet (Magorrian 2007) to compare the evolution of our 
multi-mass models against equal-mass ones. Below, we adopt A'^- 
body units G = M = a — 1. 

3.2.1 How well is the acceleration field reproduced? 

An important unsolved problem is how best to estimate the ac- 
celerations (10) given a discrete A'^-particle realization of the un- 
derlying DF /. The most sophisticated approaches to this problem 
(e.g., Dehnen (2001) and references therein) have focused on find- 
ing softening kernels that minimize the errors in the acceleration 
field given a static distribution of N equal-mass particles. In the 
present paper we do not investigate how different softening lengths 
or softening kernels affect our multi-mass models. We simply adopt 
a nested series of boxes with boundaries at \x\ — 100 x 2^* with 
i = 0, . . . , 20, each box covered by a 60"^ mesh. As one moves to 
smaller length scales the effective softening length decreases, with 
e^in = 200/60 X 2-2° ^ W^^. 

Figure 4 shows the fractional error in the radial component 
of the acceleration field returned by Grommet, in addition to the 
fractional error in enclosed mass. There is an approximately con- 
stant offset between these two quantities for equal- and multi-mass 
realization. Since our ICs here have been tailored to minimize the 
variance in the monopole component of the acceleration field, how 
important is our neglect of the higher-order multipoles? 
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Figure 4. RMS fractional deviations in acceleration (solid line), in mass 
(circle-solid line) together with its analytical value SM in dashed line; top 
panel for equal-mass model and bottom panel for multi-mass model. For 
this single acceleration calculation, we include 16 levels of refinements and 
therefore all the values should be believable outside emin ~ 10^^. 



In terms of multipole moments, the radial component of the 
acceleration is (e.g., BT87) 

Y," 



i+i r , , 



)r'''''^dr' + Ir 



J T 



where 

P!m(r)= r^dct) rdesmeYr*ie,(i,)p{r,e,4>). 

Jo Jo 
This can be rewritten as 



(31) 



(32) 



(33) 



where {Mi,n{r)) are given by 

{M,„(r)) = J fo{w')Mim{r,w')A''w' 

with projection kernels 



(34) 



Mi^{r,w') = Yr{e'A') 



--1vm(«' ) - 2; + 1 r'^+i ) 



21 + 



(35) 



where V{r) encompasses all phase-space points having radii less 
than the (real-space) radius r, and V{r) is its complement. For our 
spherical galaxy. 



{Af,„(r))= j fo{w')PIim{r,w')d^w' 



M(r), ifZ = m = 
0, otherwise. 

(36) 



The corresponding variance in 0, (7") for an A'^-body realization 
drawn from some choice of fs is 

Var {a,(r)) = 1^ ^ Var {M,„(r)) Yr{e 

Im 

where, from (17), 

Var(M,„(r)) = 1[| !MlM,^{r,w') 



2 j6 / 

d w 



(37) 

(38) 
(39) 



Similarly, the variance in the tangential component of acceler- 
ation field can be achieved by using projection kernels 

Var {aeAr)) = ^ E <ML(r)) [Yr{e, 4>)ig^^ , (40) 

Im 

where 

MUr,w')=Yr{e\(b') 



21 + 



(41) 



So, given any choice of fs, we can use the expressions above 
to calculate the contribution of the higher-order multipole moments 
to the formal eiTors in the acceleration. We find that, as we progres- 
sively include more terms, our estimate of the formal Var (ar(r)) 
approaches the actual errors observed in the A'^-body realization. 

Alternatively, we can find our optimal sampling DF fs by min- 



imizing 



E E Var(M, 



1=0 m = — l 



(42) 



truncated at say Zmax ~ 2. Notice that this new S reduces to the old 
one in eq. (18) when /max ~ 0, but otherwise includes additional 
terms with I > 0, each weighted by the monopole component I = 
0. On increasing /max from to 2, the formal Var {ar{r)) increases 
but the shape of the curve remains approximately unchanged and 
there are no noticeable differences in the resulting fs. Therefore, 
our neglect of higher-order multipole moments is justified, at least 
in the present case, provided one bears in mind that the errors in the 
full acceleration field are going to be larger by an approximately 
constant factor than what one would estimate from the monopole 
component alone. 
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Figure 5. Inner density profile of a multi-mass (top set of curves) and an 
equal-mass (lower, offset by 10 vertically) realizations of the same Hern- 
quist model evolved for 200 time units using the Grommet A^-body code. 
The ICs in each case are plotted as the heavy solid curves. The dashed 
curves show profiles after the model has been been evolved for 200 time 
units. 



3.2.2 How well are integrals of motion conserved? 

This paragraph describes the details of a full A^-body implementa- 
tion. Using both equal and multi-mass schemes, we draw 10*' par- 
ticles with radii between 10^'^ < r < 10'^. In order to suppress 
slightly deviance from symmetry (the odd terms of higher-order 
multipoles) and remove any intrinsic transient in linear momentum 
(see also McMillan & Dehnen 2005), ICs (x, v) are extended to in- 
clude the mirror distribution by reflecting each of the 10® particles 
with {x, v) -> {-X, -v). The full ICs then have iV = 2 x lO" par- 
ticles. Taking the efficiency of integration into consideration, only 
a 12-level nested series of boxes each covered by a 60"^ mesh is 
used, together with a single time-step of 2 x 10"''. Therefore, we 
expect our numerical results to be trustworthy at radii greater than 
a few times 10"''. 

Figure 5 plots the inner density profiles of both realizations 
after evolving each for 200 time-units (or 300 circular orbit periods 



at r = 0.01). The lack of particles at small radius 



10" 



the equal-mass realization means that the initial model is out of 
exactly-detailed equilibrium and causes the central density profiles 
to flatten. In contrast, the density profile in the multi-mass case is 
always much better behaved there. 

It is interesting to examine what is going on at the level of 
individual orbits. Both realizations begin with spherical symmetry 
and remain spherical, apart from the effects of Poisson noise. The 
amount of diffusion in the angular momentum J of each particle's 
orbit serves as a strong gauge of relaxation effects. This is com- 
plicated by the fact that many particles in isotropic models being 
considered here have J{t = 0) ~ 0. In such cases, even a small 
change in J(t > 0) would yield a large fractional change when 
measured in respect to its initial value. To circumvent this artifi- 
cial problem, for each particle we measure the change in angular 
momenta relative to its circular value at i = using 



AXt 



(43) 




10"' 



10"' 



10" 



10" 



10 



10' 



re/a 



Binning particles by energy and calculating the mean AX^ within 



Figure 6. Time-average diffusion rate SX^ (eq. 43) against energy labelled 
radius rg measured between t = and t = 200 for a multi-mass real- 
ization evolved using Grommet (thick curve), falcON (dashed) and for an 
equal-mass realization evolved using Grommet (thin curve). 



each energy bin gives us the time-averaged diffusion rate 5X'^{£). 
As shown in figure 6, both models suffer diffusion, but due to 
the enhancement of particle numbers and hence the smoothness 
of potential field in the central region, diffusion in the multi-mass 
scheme is suppressed by two orders of magnitude across the whole 
system. 

As a further test of the robustness of our multi-mass scheme, 
we have evolved our multi-mass ICs using the tree code FALCON 
Dehnen (2000) with a single interparticle softening radius of 10""', 
comparable to the finest mesh size used in the Grommet runs. The 
dashed curve in figure 6 plots the resulting SX^; our scheme works 
just as well for tree codes as it does for mesh codes, although the 
variable softening in Grommet does slightly decrease the amount of 
diffusion. This is not surprising, since the only difference between 
the two runs is the approximations used to estimate the accelera- 
tions (10). 

In any model with a broad spectrum of particle masses, a natu- 
ral question is what happens if heavy bodies from the outskirts visit 
the centre full of light mass elements and vice versa. To address this 
issue, we have measured the 5X^{£) of equation (43) but, instead 
of considering all particles of a given £, we compare the diffusion 
of particles on radially-biased orbits with X^ < 0.1 to those on 
nearly-circular orbits with X^ > 0.9. As shown in figure 7, there 
are no systematic differences between them. The reason for this is 
simply that particles with X ~ spend most of their time at apoc- 
entre, the apocentre radius being only a factor ~ 2 larger than the 
radius of a circular orbit of the same energy. Nevertheless, a par- 
ticle with X ~ will affect all of the more tightly bound orbits 
as it plunges through the centre of the galaxy, but our measured 
diffusion rates account for this. 



4 CONCLUSIONS 

We have presented a general multi-mass scheme to construct Monte 
Carlo realizations of collisionless galaxy models with known 
steady-state DFs /o. The scheme uses importance sampling to find 
the tailored sampling DF fa that minimizes the sum of mean-square 
uncertainties in a given set of observable quantities of the form (16). 
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Figure 7. 5X^ in multi-mass simulations for particles with < 0.1 
(dashed curve) and > 0.9 (solid). 

Although our method works for any reasonably general collision- 
less A'^-body code, we note that there are three conditions that must 
be satisfied before it can be applied successfully: 

(i) The system should be in a steady state, or close to one. 

(ii) The DF fo should be quick and cheap to evaluate, either nu- 
merically or analytically. The calculation of fs is not much more 
demanding for axisymmetric or triaxial galaxies than for spherical 
isotropic models. Finding fo for such systems is, however, non- 
trivial since one rarely has sufficient knowledge of the underlying 
potential's integrals of motion, but suitable flattened DFs do exist, 
including the standard axisymmetric two-integral f{£,Lz) models 
and also rotating triaxial models such as those used in, e.g., Berczik 
et al. (2006). An alternative way of constructing flattened multi- 
mass realizations would be to apply Holley-Bockelmann et al. 
(2002)'s adiabatic sculpting scheme to a spherical A^-body model 
constructed using our scheme. 

(iii) Finally, the utility of our multi-mass scheme depends criti- 
cally on the selection of the projection kernels Qi{w). 

Point (ii) is a well-known and longstanding problem, but the final 
condition is new. It is probably best addressed by experimenting 
with different sets of kernels, especially since it is easy to test the 
consequences of modifying them. Nevertheless, there are cases in 
which modest physical insight offers some guidance on choosing 
the Qi. Here are some examples. 

Galaxies with central massive black holes It is straight- 
forward to extend our treatment of self-consistent galaxy models to 
models containing a central black hole (hereafter MBH). By choos- 
ing kernels (as in section 3) to measure the monopole component of 
galaxy's force field and choosing fs to minimize their mean-square 
fractional uncertainty, one also achieves better spatial and mass res- 
olution within the sphere of influence of the black hole. 

Loss-cone problems The rate of supplying stars into 
MBH's loss-cone is an important ingredient in galaxy models with 
central MBHs. A thorough understanding of collisionless loss-cone 
refilling mechanisms and accurate estimates of the resulting refill- 
ing rates are particularly critical for the prediction of astrophysical 
quantities such as the timescale of binary MBH merger (Begelman 
et al. 1980; Yu 2002; Milosavljevic & Merritt 2003), the tidal dis- 
ruption rate of stars (Syer & Ulmer 1999; Magorrian & Tremaine 



1999; Wang & Merritt 2004). When using A^-body simulations to 
study such loss-cone problems, one is often interested in stars on 
low angular momentum orbits and can therefore choose kernels to 
pick out such loss-cone phase-space for detailed modelling, while 
simultaneously maintaining accurate estimates of the galaxy's ac- 
celeration field. 

Sinking satellites Kazantzidis et al. (2004) demonstrate 
the significance of using equilibrium TV-body realizations of satel- 
lite models when investigating the effect of tidal stripping of CDM 
substructure halos (satellites) orbiting inside a more massive host 
potential. Besides the shape of the background potential and the 
amount of tidal heating, the mass-loss history is very sensitive to 
the detailed density profile of the satellite itself. One can therefore 
make one step further from equal-mass realizations by designing 
kernels to pick out orbits that pass through the tidal radius, while 
again maintaining an accurate estimate of the satellite's accelera- 
tion field. 
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